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ABSTRACT: In this paper we describe spectral transformation algorithms for 
the computation of eigenvalues with positive real part of sparse nonsymmetric ma- 
trix pencils (J, L), where L is of the form ( ^ ^\ For this we define a different 



, 0^ 

extension of Mobius transforms to pencils that inhibits the effect on iterations of 
the spurious eigenvalue at infinity. These algorithms use a technique of precondi- 
tioning the initial vectors by Mobius transforms which together with shift-invert 
iterations accelerate the convergence to the desired eigenvalues. Also, we see that 
Mobius transforms can be successfully used in inhibiting the convergence to a 
known eigenvalue. Moreover, the procedure has a computational cost similar to 
power or shift-invert iterations with Mobius transforms: neither is more expensive 
than the usual shift-invert iterations with pencils. Results from tests with a con- 
crete transient stability model of an interconnected power system whose Jacobian 
matrix has order 3156 are also reported here. 

KEY WORDS: eigenvalues, stability, Mobius transforms, generalized eigenvalues 

RESUMO: Neste artigo, descrevemos algoritmos baseados em transformagoes 
espectrais para computagao de autovalores com parte real positiva de pencils de 

matrizes esparsas e nao simetricas, (J,L), em que L e da forma «). Para 

isso definimos uma extensao das transformacoes de Mobius a pencils que inibe a 
atuagao do autovalor infinito sobre as iteragoes. Esses algoritmos usam uma tecnica 
de precondicionamento dos vetores iniciais via transformadas de Mobius que junto 
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com iteragoes tipo potencia inversa com shift aceleram a convcrgencia para os 
autovalores desejados. Vemos tambcm que as transfer madas de Mobius podem 
ser usadas com sucesso no processo de inibir a convergencia para um autovalor ja 
conhecido. Alem disso, esse procedimento tem um custo computacional semelhante 
ao custo computacional de iteragoes tipo potencia ou potencia inversa com shift: 
tao caro como iteracoes tipo potencia inversa com shift aplicadas em pencAls. Sao 
tambem apresentados aqui resultados de testes com um modelo pratico para o 
problema de estabilidade transiente de um sistema de potencia interconectado, 
cuja matriz jacobiana e de ordem 3156. 

PALAVRAS-CHAVE: autovalores, estabilidade, transformagoes de Mobius, au- 
tovalores generalizados 
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1 Introduction 

The power system eletromechanical stability problem can be described by a 
nonlinear system of differential and algebraic equations 

f f{x,y) , . 

wfiere x, the state vector, contains the dynamic variables and y, the algebraic 
variables. After linearization around a system operating point (xo,yo)) i-©) 
(2^0) Z/o) such that f(xo,yo) — 0, equation (1.1) becomes 



) U3 JJ 



(1.2) 



By eliminating the vector Ay in (1.2) we obtain 

Ai = A/\x, (1.3) 

where A — Ji — J2J4^^Jz represents the system state matrix, whose eigen- 
values provide information about the singular point local stability of the 
non-linear system. The symbol A used to represent an incremental change 
from a steady-state value will be omitted from now on. 

By a classical result of ODE theory, the local stability of the system (1.1) 
can be predicted from the system (1.3). If A is diagonalizable the solution 
of (1.3) is a sum of vectors of the type e^'-^Vi, where Vi is an eigenvector 
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associated with the eigenvalue Aj. Thus, if all eigenvalues have negative real 
part, the solution decays to zero and in this case, the system is called stable. 
For an eigenvalue with positive real part, the absolute value of this expression 
increases in time and the system is unstable — these eigenvalues are called 
unstable modes. Eigenvalues with null real part give rise to oscillation, which 
never disappears. Also, eigenvalues with negative real part and non-zero 
imaginary part, but with small ratio between the real and imaginary parts, 
cause an oscillation which takes a long time to disappear — these are the 
low damped modes of the system. In the power system stability problem, 
we consider a mode A to be low damped if |i?eA| < 0.02|/mA|. In this 
paper, we search for algorithms which solve the local stability problem (1.1) 
by computing eigenvalues: we search for the eigenvalues of A with positive 
real part. The state matrices A are real, non-symmetric and dense, usually 
too large for the computation of eigenvalues by the QR method. On the 
other hand, the Jacobian matrices J are sparse and linear systems with J 
may be solved by variants of Gaussian elimination. The eigenvalue problem 
for A, 

Ax = Ax, (1.4) 



can be stated equivalently in terms of the Jacobian matrix J 
so that 



Ji J2 
J3 J4 



-^3 -hJ \yJ vO, 

or, in matrix notation, 

Jz = XLz, (1.6) 
where L = ^ ^ ) singular diagonal matrix with 1 (resp. 0) at 

diagonal entries related to state (resp. algebraic) variables and z = (^"^^ is 

an eigenvector of (J — AL). Since 2; 7^ 0, A is said to be an eigenvalue of the 
pencil (J, L) = {J — aL\a e C}. A possible approach to search for unstable 
modes is to use the shift-invert transform 



(J - \kL)zk+i = Lzk, (1.7) 

with initial shifts on the imaginary axis [16]. Although in [16] the problem is 
not described as the generalized eigenvalue problem (1.6), the system (1.7) is 
solved in order to implicitly calculate {A — XI)~^Xk — in the authors words. 
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they make use of the augmented system associated to J to shift-invert the 
state matrix A. Other methods hke subspace iteration and Arnoldi methods 
were also adapted to augmented systems and some resuhs of their apphca- 
bihty in power systems are reported in [24j . The use of the Cayley transform 
technique in order to find rightmost eigenvalues of a non-symmetric matrix 
was probably first reported in 1987 at the IEEE PICA Conference. There, 
Uchida and Nagao proposed to search for the biggest eigenvalues in absolute 
value of S" = {A + hi) {A — hl)~^, h > to detect unstable modes of the 
state matrix A of a power system [23]. There, however, the matrix- vector 
multiplication Sv was performed in a rather innefficient way. In [Ij, this op- 
eration was better implemented and the Cayley transform was extended in 
two different ways, described in detail in §2, in order to solve the power sys- 
tem stability problem, given as a generalized eigenvalue problem Jz = XLz, 
where J is not symmetric and L is diagonal with elements 1 or along the 
diagonal. The use of the Cayley transform technique to find rightmost eigen- 
values of the problem Jx = XLx, for non-symmetric J and L = ^ ^ ^ 

is also found in Computational Fluid Dynamics [Z]; an overview of this 
technique in several areas is [T7]; also, in ARPACK [13], Arnoldi iterations 
can be performed with Cayley transforms. All these references make use of 
the extension (J, L) i— )■ {J — aL)^^{J+aL), which is one of the two extensions 
of the Cayley transform analysed in [1], which turns out not to be the best 
for the problem of interest. Indeed, this extension requires a non-obvious 
strategy to control instability caused by the spurious eigenvalue at infinity of 
the generalized eigenvalue problem, in the case when L is singular. Here we 
propose to consider another extension: {J,L) i— )■ {J + aL){J — crL)^^. For 
this matrix the eigenspace associated with the eigenvalues that correspond to 
finite eigenvalues of the original problem is just the range of L, as we shall see 
in §2. Thus, the spurious eigenvalue is handled by keeping iterations in this 
space. In the power system stability problem, the range of L is the set of vec- 
tors with coordinates related to algebraic variables equal to zero. In the case 

L=^Q Q^isa matrix of order n + m, where M is a x n positive definite 

matrix, as in [6], the space is just the set of vectors with the last m coordi- 
nates equal to zero. In §2, we introduce Mobius transforms (a generalization 
of the Cayley transforms) for the generalized non-symmetric eigenvalue prob- 
lem. Mobius transforms can be used to precondition random initial vectors 
as well as to inhibit the convergence to eigenvalues already found (much like 



4 



a deflation technique), at a computational cost not more expensive than a 
shift-invert iteration with pencils. The application of polynomial fllters to 
vectors in the computation of eigenvalues of sparse non-symmetric matrices 
has been the subject of several papers [21], [18], [5]. We will see that Mobius 
transforms can be seen as the action of a rational function fllter which gives 
infinite and zero weights to two arbitrary points in the complex plane. These 
techniques are presented in §3 together with two algorithms based on them. 
Finally, in §4 we apply four methods to the computation of the unstable 
modes of a pencil (J, L), where J is a sparse matrix of order 3156 and L is 
diagonal with elements 1 or on the diagonal, of rank 790: the algorithms 
we suggest, an implementation of the Arnoldi method (ARPACK) and an 
implementation of the Arnoldi method with acceleration by Chebyshev poly- 
nomials (ARNCHEB), both applied on the extension of a Cayley transform 
introduced here and a subspace iteration applied on the pencil {J,L). 

2 The Generalized Eigenvalue Problem and 
Mobius Transforms 

The Mobius transforms are complex functions 

Ck,a,p ■ C U {OO} > C U {OO} 

s ^ k ^ 

s — a 

where a + (5 ^ 0. They are conformal mappings which map lines and 
circumferences to lines or circumferences. When k = 1 and /3 = a, with 
Re a 7^ 0, these functions are the so called Cayley transforms, which map 
the semiplane {z\Rez > 0} to {z;\z\ > 1} ({2; < 1}) if Rea > 
{Re a < 0). Mobius transforms can be defined in the space of square ma- 
trices in a analogous way: given a square matrix A and a not in X{A) (the 
spectrum of A) then k{A + PI){A — al)~^ is a matrix whose spectrum is 
{k{X + (3){X — a)~^; A G X{A)}. However, the extension of these transforms 
to pencils (J, L) can be done in two ways. The goal of this section is to 
present the advantages of one extension over the other in the application 
to the power system stability problem. Also, we will see that multiplica- 
tion of a Mobius transform against a vector requires no more computations 
as solving the equation (J — aL)z = w, for some (easily computed) scalar 
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a. We begin with a lemma that insures that under generic conditions the 
pencil (J, L) has exactly n eigenvalues, where n is the rank of L. Since 
by the Singular Value Decomposition there are unitary matrices U, V such 
that U{J — aL)V^ = J — aJ2, S = diag{di, dn, 0, 0) for appropriate real 
numbers [9], we may suppose without loss that L = diag{di, dn, 0, 0). 

Lemma 21 Let J ~ {^j ^ ^ ~ o)' ^^^'^^ D is a non- 

singular diagonal matrix. Then, if J4 is nonsingular, (J — aL) is singular if 
and only if a is an eigenvalue of D^^A, where A = Ji — JiJ^^Jz- 



Proof: If J4 is nonsingular, A = Ji — J2-T^'^3 can be defined. Now, 

=^ (A — aD) X = and y = —J^^J^ x. 
Therefore, (J — aL) is singular <^==^ a is an eigenvalue of D~^A. 



Ji - aD J2\ fx\ _ fO 
J3 Jj[y)~[0 



From now on, the non-sigularity of J4 will be assumed. Let a G C, such 
that (J — aL) is nonsingular. For k, P & C, k ^ and (3 7^ —a, let 

Ck,a,i3 = k{J + (3L){J - aL)-^ and Dk,a,f3 = k{J - aL)~\j + (3L). 

Notice that these two matrices have the same spectrum. Moreover, they are 
both extensions of Mobius transforms applied to the pencil {J,L) and the 
spectrum of (J, L) is related to the spectrum of the two extensions by the 

s + ^ 

relation s k . However, the eigenspaces of both extensions are not 

s — a 

the same in general, according to the following propositions whose proofs are 
left to the reader. 

Proposition 22 (a) z ^ is an eigenvector of Ck,a.i3 associated with fi, 

H ^ k, if and only if Jw = XLw, where w = (J — aL)^^ z and X = ^ ~^ ^ — _ 

fj, — k 

(b) z ^ is an eigenvector of Ck^a,i5 associated with the eigenvalue k if 
and only if L{J — aL)^^z = 0. 

Proposition 23 (a) z ^ is an eigenvector of D^^a^p associated with fi, 

fi ^ k, if and only if Jz = XLz, where X = ^ . 

fi — k 

(b) z ^ is an eigenvector of Di^^a,i3 associated with the eigenvalue k if 
and only if Lz = 0. 
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Thus, the finite eigenvalues of the pencil {J,L) correspond to the eigen- 
values of Ck,a,ii or Dk^a,i3 that are different from k. These eigenspaces are 
described in the following proposition. 

Proposition 24 (a) The eigenspace of 0^^^,^ associated with the eigen- 
values different from k is the range of L. 

(b) The eigenspace of Dk,a.iB associated with the eigenvalues different from 
k is the range of (J — aL)^^L. 

Proof: Let ^ be a vector such that C^f^^^^z = kz. Thus (J^ — (5L^)z = 
{J^ — aL^)z and, since /3 a, L^z — 0. Since the eigenspace of Ck,a,i3 as- 
sociated with the eigenvalues different from k is orthogonal to the eigenspace 
of C^^ p associated with A;, that eigenspace is the range of L. 

The proof of the second part of the proposition is analogous. ■ 



Corollciry 25 Let J — ( j ) a matrix of order n + m, where J4 is 



nonsingular matrix. Then the eigenspace of Ck,a,f5 associated with the eigen- 
values different from k is the space of vectors whose last m coordinates are 
zero. 

Since we are interested in calculating eigenvalues of the pencil (J, L) from 
Mobius transforms, the eigenvalue k of the transform must be treated with 
special care. The corollary above states that in the case of the power system 
stability problem, for instance, the invariant subspace V corresponding to the 
eigenvalues different from k is the set of vectors v that have null coordinates 
in the positions related to algebraic variables. The iterates of our approx- 
imate eigenvectors ought to stay in this subspace, and if the computation 
of Ck^a,i3V leaves V because of errors due to finite precision arithmetic, we 
simply project the results back to V by zeroing the appropriate coordinates. 
The analogous iteration with the Mobius extension Dk^a,p is not subject to 
such an easy stabilization procedure, and in this case approximate eigenvec- 
tors will frequently converge to the eigenspace associated to k, which is of no 
real interest for the pencil eigenvalue problem. Hence, we will only consider 
here iterations making use of the extension Ck,a,i3- 

Now, let (T e C, with i^ecr > 0, such that (J — aL) is nonsingular, and 
consider 




where M is a n x n 



-1 



(2.1) 
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Let 7^ 1. Then, simple calculations obtain 

{C^-liI) = {I - + 2Rea L{J - a L)-^ 

-L jJi jJj J- 1^ -L 

Thus, the computational cost of applying any of the three matrices in the 
left-hand side to a vector is equivalent to solving a system (J — aL)w = z. 
Similar results hold for both extensions of Mobius transforms. 

The spectra of the pencil (J, L) and of C^. are related to each other by 
the bijective function 

c^:CU{oo} — )■ CU{oo} 
s + W 

s 1-^ 

s — a 

This function maps complex numbers with negative real part onto the unitary 
circle, pure imaginary ones onto the unitary circumference and those with 
positive real part onto numbers outside the unitary circle. The search for 
eigenvalues of largest absolute value of the extension of the Cayley transform 
thus obtains the unstable modes of the original pencil. Unfortunately, this 
transformation clusters some eigenvalues very close together, thus affecting 
adversely the rate of convergence. Therefore, we need techniques to accelerate 
the convergence to the desired eigenvalues. This is the goal of the next 
section. 



3 A Class of Spectral Algorithms 

Two ways of extending Mobius transforms to the pencil (J, L) were described 
in the previous section. When L is a matrix of order n + m and of the 

type ^ Q Q ^ ; where M is an n x n nonsingular matrix, the eigenvectors 

of the extension Ck,a,i3 associated with eigenvalues different from k are the 
vectors which have zero entries in the last m coordinates. Therefore, in this 
case, the attraction of the eigenvalue /c, which corresponds to the infinite 
eigenvalue of the pencil, can be easily avoided. Now, the power system 



8 



stability problem, where M is the n x n identity matrix, has additional 
features that lead us to explore Mobius techniques. Usually, this problem 
has several negative real eigenvalues with large modulii, which are mapped 
to eigenvalues close to k: the clustering of eigenvalues substantially reduces 
the speed of convergence of the power method [9] . In this section we introduce 
a class of algorithms which use Mobius transforms to precondition vectors, in 
order to yield more convenient initial vectors for a search process with shift- 
invert Mobius transforms iterations inside the unit circle. Also, we introduce 
a way of inhibiting known eigenvalues in the iteration by yet another use of 
Mobius transforms. 

3.1 Preconditioning and Shifts 

Let a E C, Rea > 0. As seen in §2, the eigenvalues of {J,L) with positive 
real part correspond to the eigenvalues of located inside the unit circle. 
In order to achieve larger components of eigenvectors of associated with 
eigenvalues of modulus less than one in the search vectors, we start with 
random vectors and apply C^- to them a few times, obtaining the so called 
preconditioned vectors. We then use these vectors in a search process for 
eigenvalues inside the unit circle. How should one choose a ? The reality of 
the pencil (J, L) implies that eigenvalues come in conjugate pairs, which is 
still true for the matrix C„ if a is taken to be real: our search for eigenvalues 
is then reduced to, say, the upper half-disk. Also, in this case, a simple 
differentiation shows that 

a = \X\ 

maximizes |ca-(A)| = where X = a + bi, b ^ — thus, the choice 

(7 = |A| takes eigenvalues of the pencil of absolute value |A| to eigenvalues 
of Ccr of largest possible absolute value. We then take shift-invert iterations 
with as the matrix to be shifted. Since the eigenvalues of are the 
inverse of the ones of Co-, the dominant eigenvectors of C^- are associated 
with the eigenvalues of which are either inside the unit disk or near the 
unit circle. Based on these remarks we introduce the first algorithm. 

• Algorithm I 

Step 1 Begin with r orthogonal vectors belonging to Af{L)-^. 
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Step 2 Multiply the vectors by p times, normalizing them after each 
multiplication (e.g., keep them with sup norm equal to one). Let 
f j, ? = 1, r, be the resulting vectors. 

Step 3 Take initial shifts /c = 1, s, in a circumference of radius e, 
< e < 1. 

Step 4 Let wf^ = Vi, i = 1, ...,r. 
For k = 1, ...,s 

fori = 1,2,... 

= u^p /af \ where is the coordinate of of 
maximum absolute value; 

ifj=t,2t,... 

= fi^^ H ^ , where q is such that 

\\wy^ — wi^'^^W^ = min ||i(7p^ — Wr''~^^\' 



■■q "^q 



loo iV-"-^ W^i "^i Moo- 



Preconditioning is performed in Step 2 — different choices of p are pre- 
sented in the experiments of §4. The shifts in Step 3 are taken in a circle 
centered at the origin. The convergence of shift-invert iterations to an eigen- 
value depends enormously on the choice of shift: the aim of this shifting 
strategy is to cover the unit circle. It is at this point that the choice of a 
real parameter a entitles us to divide by two the search for eigenvalues by 
taking into account the reality of the original pencil. Every t iterations, the 
shift /i? is updated by the formula above ^25j . When the variation of one of 
the vectors between the previous and the current iterations is smaller than 
a fixed tolerance, the shift is updated and the resulting value is taken to be 
an eigenvalue. 



3.2 Inhibiting Convergence of Eigenvalues 

Mobius transforms could be used to inhibit the convergence to an eigenvalue 
^ of if in Algorithm I iteration vectors were multiplied by {C^^ — ^I). 
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Notice that after t iterations with the matrix (C^ ^ —^I) (C^ ^ (again 
a Mobius transform) the updated shift should be 

4^ )-i ' 

where o;^-'^ is defined in a similar way. The computational cost is again no 
more expensive than a standard shift-invert step for the generalized eigen- 
problem because 

The algorithm below uses this deflation-type strategy. 
• Algorithm II 

Use the Algorithm I to obtain preconditioned vectors wf^ — Vi, i — 
1, r, and initial shifts k — 1, s. 

For A; = 1 

Follow Algorithm I up to convergence to an eigenvalue ^. 
For A; = 2, 

for j = 1,2,... 

= u^p /a^P ^ where is the coordinate of of max- 
imum absolute value; 

if J =i,2i,... 

l^k TV ' where q is such that 

a^q - 1 

I \wl^^ — W^J~^^ I loo = Klin I — w\^~^^ I loo- 
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Figure 1 : Sparse Pattern of the Test Matrix 

4 Tests and Comparisons 

We now present some results of experiments made with the algorithms de- 
scribed in the previous section, together with tests performed with a subspace 
iteration method ([H], [22]) and an Arnoldi method ([2], [3] [H], [19]). The 
test pencil (J, L) is taken from a transient stability model of the South- 
Southeast interconnected Brazilian power system: J is the Jacobian matrix 
at an operating point and is of order 3156 and L is a diagonal matrix with 
elements 1 or 0, corresponding respectively to state and algebraic variables, 
with rank equal to 790. Only 0.14% of the elements of J are nonzero and its 
sparseness pattern is given in Figure 1. This pencil has exactly four eigen- 
values with positive real part: 0.1814±i4.8323, 0.0233 and 0.0004. The first 
two correspond to genuine unstable modes of the system. The remaining 
two are related to two redundant states of the system [12] : they would have 
been zero if there were no roundoff errors in the generation of the Jacobian 
matrix. 

The same pencil has been used in [1], where the authors report re- 
sults of a parallelization of the lopsided simultaneous iteration algorithm 
^22j . Their strategy to find unstable modes is to perform shift-invert sub- 
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converged value 


iter 


prod 




converged value 


Iter 


prod 


-0.1164+3.2018i 


3 


96 




-0.0925+3.98271 


2 


64 


-0.0925+3.9827i 


6 


180 




-0.1164+3.20181 


7 


204 


-0.9573+2. 1594i 


8 


228 




0.1814+4.83231 


7 


204 


-0.4803+1. 7632i 


8 


228 




-0.5911+4.69351 


8 


224 



Table 1: Simultaneous shift-Invert method: a = + 3i and a = + Ai 

space Iterations with Initial shifts given on the Imaginary axis. We Imple- 
mented here a sequential version of this algorithm. As Initial vectors we 

took Z^^^ = LZ^^^ = i Q j, where X^^^ Is a matrix with (up to) eight 

first column vectors of the 790 x 790 Fourier matrix F{j, k) = e-^ '^'^. After 
each four shift-invert iterations followed by a normalization of the vectors a 
Rayleigh-Rltz acceleration was done. That is, we first calculated a spectral 
decomposition oi B = Q-^H ([13j, [20]), where G = {LZ^^'^YiLZ^^''^) and 
H = {LZ^^^^)^W, with W = L{J - aL)-^LZ^^^\ such that the Rltz values 
were ordered according to decreasing absolute values. Then W was multi- 
plied by the matrix of Rltz vectors and the resulting column vectors, after 
normalization, undertook another cycle of four Iterations. Convergence was 
achieved when the oo — norm of the difference between corresponding vectors 
of LZ(^^) and LZ^^'^+i) was less than a tolerance, taken to be 10 ^ . A defla- 
tion technique was also Implemented for these tests: aside from the choice of 
initial vectors and the deflation technique, this is the algorithm used in [1]. 
Some results obtained with these Iterations on SUN SPARC workstations 
are displayed in Table 1, where iter and prod mean respectively the number 
of Rayleigh-Rltz accelerations and the number of matrix-vector multiplica- 
tions being performed. Notice that matrix-vector multiplication consists of 
backward and forward substitutions after computing an CLi decomposition 
of (J — aL), which is done only once for each a. 

ARNCHEB performs an Incomplete Arnoldl method combined with an 
acceleration technique using Chebyshev polynomials |3]- We have used its 
Reverse Communication interface to calculate matrix-vector multiplications 

by Ca = I + 2Rea L{J — oL)^^. As seen before, since L = ( [| V we 



,0 0. 

avoid the spurious eigenvalue 1 with this Cayley transform by considering 
only the first 790 coordinates of the vectors, that is, by solving systems 
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converged value 


0(V) 




converged value 


0(V) 


0.1814+4.83231 


10-10 




0.1814+4.83231 


10-11 


0.1814-4.83231 


10-10 




0.1814-4.83231 


10-11 


0.0233+0.00001 


10-08 




0.0233+0.00001 


10-08 


0.0004+0.00001 


10-07 




0.0004+0.00001 


10-08 


iter 


prod 




iter 


l)rod 


12 


856 




31 


2107 


Table 2: ARNCHEB: a 


= 4.0 and a = 6.0 


converged value 


0(V) 




converged value 


0(V) 


0.1814+4.83231 


10-13 




0.1814+4.83231 


10-14 


0.1814-4.83231 


10-13 




0.1814-4.83231 


10-14 


0.0233+0.00001 


10-14 




0.0233+O.OOOOi 


10-14 


0.0004+0.00001 


10-13 




0.0004+0.00001 


10-13 


Iter 


prod 




Iter 


prod 


301 


3349 




301 


3679 



Table 3: ARPACK: a = 4.0 and a = 6.0 

(J — aV)w — Lz and taking only the first 790 coordinates of In the tests, 
the dimension of the Krylov space was taken to be 54 and the number of 
requested eigenvalues, 4. Some results are In Table 2, where 0{W) Is the 
order of the residual \ \{Ca — \I)v\\2i with 11^112 = 1, and iter Is the number 
of Arnoldl steps. 

The same sort of experiment was carried out with ARPACK. The program 
dndrvl.f was rewritten to Include the same routines which solved the systems 
In the tests with ARNCHEB. Here the dimension of the Krylov space was 
taken to be 20 and also 4 eigenvalues were requested to converge. Some of 
the results are shown In Table 3. 

The advantage of these methods Is that they require only one factorization 
of the pencil In upper and lower triangular factors. The disadvantage Is the 
presence of parameters which need to be adjusted, like the dimension of the 
Krylov space and the convergence criterion. ARNCHEB worked well when 
this dimension was large compared to the number of desired eigenvalues 
(ten times, e.g.). The tolerance used was l.Od-11. The opposite ocurred 
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converged value 


Iter 


0(V) 




converged value 


Iter 


0(V) 


0.0004+O.OOOOi 


6(2) 


10-10 




0.0004+O.OOOOi 


6(2) 


10- 


-10 


-0.6223+0.96491 


11 (3) 


10-10 




0.0233+0.00001 


10 (3) 


10- 


-11 


-0.4803+1.76321 


8(2) 


10-07 




0.1814+4.83231 


10 (3) 


10- 


-08 


-0.0925+3.98271 


6(2) 


10-11 




-0.0925+3.98271 


6(2) 


10- 


-09 


-0.1351+6.89741 


8(2) 


10-07 




-0.1351+6.89741 


9(3) 


10- 


-11 


-1.5684+12.5931 


7(2) 


10-08 




-0.1144+10.6171 


8(2) 


10- 


-06 


Table 4: Algorithm I ■ 


- a — 


4.0: p = and p = 


= 40 







with ARPACK: It worked well when the dimension of the Krylov space was 
between four and six times the number of desired eigenvalues. The tolerance 
employed was O.dO: when It was changed to l.Od-16, the convergence for the 
two positive real eigenvalues was not achieved. 

The experiments with algorithms I and II were carried out by taking only 
the first four column vectors of the Fourier matrix of order 790 as Initial 
vectors. We chose Initial shifts on the upper half of the unit circle, /Xfe = 
e~^, 0<A;<6(A; = 6 yields (C^i — /), which Is singular). For each 
Initial shift, the process stops when, for some /c, \x{^ — x^i^^\ < tol. The 
tolerance tol was taken as 10-^. The number t of Iterations performed before 
shift update was 4. The following tables contain the results of these tests 
performed In SUN SPARC workstations. The first column Indicates the 
eigenvalue reached by the Iteration. Also, iter Is the number of shift-Invert 
Iterations, the number of JCU factorizations appears between parenthesis, and 
0(V) Is the order of the residual 

||(C;1-/X/)X||2, |k||2 = l. 

Preconditioning Is performed by multiplying Initial vectors p times by C^, 
followed by normalization: p — means no preconditioning. If p ^ 0, 
multiplication by requires an additional CU decomposition. In order to 
know how many matrix-vector multiplications were performed In the tests, 
multiply iter by 4 and add 4p. 

For the algorithm II, a was chosen to be 4.8334, the modulus of the 
unstable eigenvalues 0.1814 ± 4.8323i. From the previous section, this Is 
the real value for a that maximizes the absolute value of the corresponding 
eigenvalues of the Cayley transform Ca- Thus, preconditioning of the Initial 
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A 


Iter 


0(V) 




converged value 


Iter 


0(V) 


0.0004+O.OOOOi 


6(2) 


10-10 




0004+n nnnoi 


6 (2) 


10-10 

X w 


-0.4803+1. 7632i 


7(2) 


10-07 




0233+0 00001 


10 (3) 


10-08 


-0.1164+3.2018i 


8(2) 


10-07 




n 1814-1-4 8323i 


10 ('31 


10-08 

X w 


-0.1764+6.1231i 


7(2) 


10-08 




-0 1764-1-6 1931i 




-10-07 


-0.1144+10.617i 


5(2) 


10-10 




-D 1 1 44-1-1 n 61 7i 

W . X X tit: X W . W X ( 1 


4 fTl 


10-06 

X w 


-1.5684+12.593i 


11 (3) 


10-09 




_1 01 Q^-l-n OOOOi 


18 f5l 


10-09 


Table 5: Algorithm I ■ 


- a — 


6.0: p = and p = 


= 40 




A 


Iter 


0(V) 




converged value 


Iter 


0{V) 


0.0004+O.OOOOi 


6(2) 


10-10 




0.0004+0.00001 


6(2) 


10-10 


0.0233+O.OOOOi 


10 (3) 


10-09 




-0.1764+6.12311 


12(3) 


10-07 


0.1814+4.83231 


12 (3) 


10-06 




0.0233+0.00001 


10 (3) 


10-09 


0.1814+4.83231 


6(2) 


10-08 




0.1814+4.83231 


6(2) 


10-07 


-0.1144+10.6171 


9(3) 


10-12 




-0.1144+10.6171 


8(2) 


10-08 


-200.38+0.00001 


15(4) 


10-09 




-200.38+0.00001 


11(3) 


10-12 



Table 6: Algorithm I and II - cr = 4.8334, p = 40 

vectors with this a should make these eigenvalues easier to detect. Indeed, 
In Table 6 we can see that one of them was Identified from two consecutive 
Initial shifts. In the same table are listed the results when the procedure 
of Inhibiting the last found eigenvalue was applied (for the first Initial shift, 
there Is nothing to Inhibit). In Table 7 a better performance of the Inhibiting 
procedure can be seen: the convergence to the two unstable eigenvalues was 
achieved. 

The jCU decomposition of (J — aL) Is about 150 times slower than the 
resolution of the corresponding systems on a SUN SPARC 4 workstation 
(typical runtimes were 7.05078 and 4.29688e-02, respectively). Thus, Algo- 
rithms I and II had a performance comparable to the others In respect to 
time and accuracy. For Instance, from Table 4, we see that by precondi- 
tioning the vectors we obtained three unstable modes and three stable ones 
(two of these are low damped), after 196+160 matrix- vector multiplications 
(Indeed, backward and forward substitutions) and 15+1 factorizations. 
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A 


iter 


0(V) 




converged value 


iter 


0(V) 


0.0004+O.OOOOi 


6(2) 


10-1° 




0.0004+O.OOOOi 


6(3) 


10-10 


0.0233+O.OOOOi 


8(2) 


10-06 




0.1814+4.8323i 


13 (5) 


10-09 


0.1814+4.8323i 


7(2) 


10-06 




0.1814-4.8323i 


9(4) 


10-08 


0.1814+4.8323i 


6(2) 


10-10 




0.1814+4.8323i 


6(3) 


10-09 


-0.1144+10.617i 


13(4) 


10-12 




-0.1144+10.617i 


7(3) 


10-07 


-0.1144+10.617i 


8(2) 


10-07 




0.1814-4.8323i 


10 (4) 


10-07 



Table 7: Algorithm I and II - cr = 4.8334, p = 80 
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